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Abstract. One of the most used approaches in simulating materials is the tight-binding approximation. When using this 
method in a material simulation, it is necessary to compute the eigenvalues and eigenvectors of the Hamiltonian describing the 
system. In general, the system possesses few explicit symmetries. Due to them, the problem has many degenerate eigenvalues. 
The ambiguity in choosing a orthonormal basis of the invariant subspaces, associated with degenerate eigenvalues, will result 
in eigenvectors which are not invariant under the action of the symmetry operators in matrix form. A meaningful computation 
of the eigenvectors needs to take those symmetries into account. A natural choice is a set of eigenvectors, which simultaneously 
diagonalizes the Hamiltonian and the symmetry matrices. This is possible because all the matrices commute with each other. 
The simultaneous eigenvectors and the corresponding eigenvalues will be in a parametrized form in terms of the lattice 
momentum components. This functional dependence of the eigenvalues is the dispersion relation and describes the band 
structure of a material. Therefore it is important to find this functional dependence in any numerical computation related to 
material properties. 
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INTRODUCTION 

Tight-binding (TB) is a method used to investigate the electronic structure of a large class of solid materials fUj. 
When used in conjunction with numerical simulations, this method introduces several simplifications that reduce the 
complexity of the description of the material. Every solid material is constituted of atomic nuclei that identify a 
lattice and are the source of potential energy. On other hand, the nucleus-nucleus and electron-electron interactions 
are neglected. The TB model assumes that the electrons are tightly bound to their corresponding nuclei, implying that 
their wave functions are localized. Furthermore, atoms interact weakly only through their valence electrons. 

Since the electrons are moving independently, the Hamiltonian H of the system is given as a sum of the kinetic 
energies of the electrons [p^ /2m) and the potentials due to the nuclei Y.i^{fi^Rn), with Rn being the positions of the 
nuclei in three-dimensional space r,. Thus: 



H = yH, = y\!^ + yv{n-R„)\^ (d 




where Ng denotes the number of considered electrons and A^j the number of lattice sites of the crystal IJl |3] . To find the 
eigenstates of this system, a linear combination of the atomic orbitals (LCAO) Lj V,^,- is used as an ansatz. The atomic 
orbitals 0, are the eigenstates of the Hamiltonian for an isolated atom and v,- the coefficients to be computed with the 
constraint that | | = 1. Since the overlap of the atomic orbitals of neighboring atoms is assumed to be small, they 

are treated as orthonormal, i.e. their inner products are (0,, 0,) = 5,^. Using this property and the LCAO as an ansatz, 
one obtains the following eigenproblem: 

Hv„^e„v„ with n = 0,1, . . . ,K — 1 , (2) 

where H £ C^^^ is the Hamiltonian in the basis of the atomic orbitals. The quantity v„ S is a vector of coefficients 
V, of the LCAO and the eigenvalue e„ e M is the associated energy level. Note that the Hamiltonian is hermitian and 
therefore the eigenvalues are real. 

The entries of the Hamiltonian are given by 



Hke = Skecck + jike, 



(3) 



where is the resuh of an overlap integral between neighboring electronic orbitals and the underlying lattice 
potential OJ. In the simplest case of equal atoms and only nearest neighbor interaction, the expressions in Eq. ([3]) 
simplify to jSjtc — — f [Skj^i + dk+ij) and aic = a, where a and f are constants |3|. Since f represents the interaction 
between neighboring atoms, it is often called the hopping term. 

The quantum mechanical problem of finding the electron wave function is therefore reduced to the solution of a 
finite dimensional eigenproblem. Having computed the eigenvalues and eigenvectors, we aim at expressing them in 
terms of the lattice momentum components k = (fci ,^2,^3)- Eventually, the whole set of eigenvalues can be seen as a 
function of k, called the dispersion relation. This is an important relation from which we can determine a large set of 
physical properties of a material ISO. Therefore determining this relation numerically is our final goal. 



A 2-DIMENSIONAL EXAMPLE 



In this section we construct a simple example. While it can be solved analytically, we show that it can also be 
accurately solved numerically. Consider a two-dimensional rectangular lattice of equal atoms, as shown in Figure 
[T](left). The dark colored atoms constitute our N-hy-N lattice structure and the brighter atoms represent the use of 
periodic boundary conditions. Each atom in the structure interacts with its four nearest neighbors. The interaction is 
given by the hopping term t as discussed above. 




Figure 1. Left: Two-dimensional rectangular lattice of equal atoms and nearest neighbor interaction. In this figure the mesh has 
the size N = 4; right: Spectrum of the Hamiltonian with a = 1.0 and t = 0.2 for A' — 8. 



The Hamiltonian H e 



of the system has the form 
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(4) 



The matrix C <E W^^'^ is circulant. It is equivalent to the Hamiltonian for the one-dimensional lattice of identical 
atoms with periodic boundary conditions and nearest neighbor interactions. It has the same structure as H in Eq. Q 
with C and D replaced by the scalars a and — f, respectively. The matrix D E M.^^^ is diagonal with all elements equal 
to -f . 



The eigenvalues and eigenvectors of H can be expressed in closed form (see ||4]|5|). Since H is block-circulant with 
circulant symmetric blocks, its normalized eigenvectors v„ e are 



Vn{r,s) 



N 



PrWs 
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4^ 



, Pr = exp ( i^r \ , t,s = exp ^'^i ) , (5) 



where r,s = 0,1,. ..,A^ — 1. The parameters p, and ^s are the A^-th roots of unity. The eigenvectors v„ form an 
orthonormal basis ID. The corresponding eigenvalues are 



e„(r,i) = a — 2fcos ( "^'"1 — 2fcos I — i 



(6) 



All eigenpairs (e„,v„) are parametrized by the quantities r and s. The index n of the pair can be defined as any 
bijective function n — f{r,s). 

When defining kj^ and ky as k^ :— 2nr/N and ky :— 2ns/N, respectively, Eq. ( [6| descr ibes the dispersion relation 

e{kx,ky). This relation yields all the allowed energies for possible momenta k = yl^~+k^. . Because of the periodicity 

of the crystal, both energy and momentum are quantized ||2l[3|. As N ^ °°, the dispersion relation reveals the band 
structure of the crystal. Therefore it is important to identify this relation. In more complicated cases for which no 
analytical solution is available, it is important to compute the dispersion relation through a numerical procedure. 

In general a numerical computation does not result in eigenvectors in the form of Eq. (j5]l. For every m-degenerate 
eigenvalue, the matrix H only defines an m-dimensional invariant subspace. For instance, from Eq. (|6]l it can be seen 
that for even only the largest and smallest eigenvalues are distinct. All the other eigenvalues have at least multiplicity 
four Figure[T](right) illustrates this behavior showing the spectrum of a Hamiltonian for = 8. The degeneracy creates 
an ambiguity in choosing a basis for the associated invariant subspaces. 

The best a general algorithm can do given only the matrix H, is to compute an arbitrary orthonormal basis. 
Such a solution would not have a parametric expression in terms of the momentum because it does not respect the 
symmetry of the lattice. In physical terms, it means that the set of eigenvectors is not invariant under the action of the 
symmetries expressed in matrix form. An obvious solution is to find a vector basis that simultaneously diagonalizes 
the Hamiltonian H and the symmetries Si. This is indeed possible, because H commutes with the symmetric^ 12, 
generating a closed algebra under multiplicatiorj^ 



(7) 



Our example has two translation symmetries, one along the x-axis Sx = In ^ Cp, and another along the y-axis 



Sy = Cp®lN. The symbol ® denotes the Kronecker or tensor product. In € 



TiNxN 



is the identity matrix, and Cp G 



liNxN 



is a circulant matrix with Cp = circ(0, ... ,0, 1). Since both matrices are simply permutation matrices, their inverses 
are Sj'^ = Sj , which correspond to translations along the negative direction of the symmetry axis. 

Unfortunately, neither nor Sy have distinct eigenvalues, in fact every eigenvalue has multiplicity A^, and we can 
only determine invariant subspaces. In order to find the simultaneous eigenvectors of all the matrices, we can create 
linear combinations that are part of the algebra. Using enough combinations will allow us to identify uniquely (up to 
a phase) the simultaneous eigenvectors in their parametric form. 

For example, we can compute the eigenvectors of H{Sx — Sy) and Sx{II ~ Sy), and select from their eigenvectors 
only those that are simultaneous eigenvectors of all the matrices. In this way, we generate a set of A^^ simultaneous 
eigenvectors. They coincide with the analytical solutions of Eq. Q up to a phase factor In our example we can 
normahze the eigenvectors imposing the first element to be real. The resulting basis is uniquely defined and satisfies 
the symmetries of the problem. 



' This is a well-known result from Hamiltonian dynamics: the time-dependence of a generic operator A is described by the Eq. ^ <x [H,A\. If A 
represents a conserved symmetry, its derivative with respect to time is automatically null from which the thesis follow. 
^ The Jacobi identity is readily verified. 



To show a concrete numerical example, we now look at the results for a 25-by-25 lattice with a = 1 .0 and t — 0.2: 
using MATLAB® for the computation, and denoting computed quantities with a hat, the maximum error of the 
computed eigenvector entries compared to the analytical solution for real and imaginary part is about 3.5 • lO^''*; the 
maximum residual max/ || i/v, — gjV,- ||2~ 1.9- lO^'^''; the orthogonality of the computed eigenvectors is max;j |v*Vy — 
5ij \ « 3.8 • 10^'^; finally the maximum error in the eigenvalues is max,- |9{(vJ'i/v,) — e, | « 1.1 • 10^'^. 

Figure|2](left) shows the dispersion relation computed numerically, while on the right we present the error compared 
to the analytical solution given by Eq. (|6]). 




Figure 2. Left: The dispersion relation e{kx, fcy) computed numerically; right: The error compared to the analytical solution given 
by Eq. 



SUMMARY AND CONCLUSION 

The computation of eigenvalues and eigenvectors of an Hamiltonian describing a quantum mechanical system can lead 
to eigenvectors that are not satisfying certain physical requirements. 

Through a simple model of a solid material, we show that symmetries introduce degeneracies in the spectrum of 
the Hamiltonian. These degeneracies lead to an ambiguity in computing a basis for the invariant subspaces associated 
with the corresponding eigenvalues. A basis, if not chosen carefully, in general will lead to eigenstates which are not 
taking the symmetries of the problem into account. In order to generate a satisfactory basis the eigenvectors must 
simultaneously diagonalize the Hamiltonian and the symmetry operators. The eigenvectors that fulfill such a condition 
compose a complete orthonormal eigenvector basis that is uniquely defined. Finding this basis is the first step in 
computing the dispersion relation of the material under investigation. 

We explored the concept of computing numerically the dispersion relation in a simple model having analytical 
solutions. Our final goal is to apply this approach to the investigation of irregular materials, where analytical solutions 
are not known. 
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